#ifndef  QUADRATICBSPLINES_H
#define  QUADRATICBSPLINES_H

#include "Spline.h"
#include "B_Basis_function.h"
//二阶样条

class QuadraticBSplines:public BSplines{
    public:
    std::vector<double>knots;
    std::vector<double>f;
    std::vector<double>coefficients;
    double f_1;
    double f_n;

    public:
    QuadraticBSplines(std::vector<double>knots,std::vector<double>f,double f_1,double f_n);
    void CalculateCoefficients() ;
    double B_value(double x_value);
    };

    QuadraticBSplines::QuadraticBSplines(std::vector<double>knots,std::vector<double>f,double f_1,double f_n):BSplines(knots,f),knots(knots),f(f),f_1(f_1),f_n(f_n){        
        
        coefficients.resize(knots.size()+1); 
        CalculateCoefficients();
    }
void QuadraticBSplines::CalculateCoefficients(){
    std::cout<<"qudratic"<<endl;
    int n=knots.size();
    Eigen::MatrixXd A=Eigen::MatrixXd::Zero(n+2,n+2);
    Eigen::VectorXd b=Eigen::VectorXd::Zero(n+2);
    //三次样条，每个区间上都受到三个函数的影响
    for(int i = 1; i <= n; i++){
        for(int j = i-1; j <= i+1; j++){
            B_Basis_function B(2, j-1, knots);
            A(i, j) = B.value(knots[i-1]);
        }
        b[i] = f[i-1];
    }
    
    //边界条件 
    B_Basis_function B__1(2,-1,knots);
    B_Basis_function B_0(2,0,knots);
    B_Basis_function B_1(2,1,knots);

    B_Basis_function B_n__2(2,n-2,knots);
    B_Basis_function B_n__1(2,n-1,knots);
    B_Basis_function B_n(2,n,knots);



    A(0,0)=B__1.value(-11/2.0);
    A(0,1)=B_0.value(-11/2.0);
    A(0,2)=B_1.value(-11/2.0);
    A(n+1,n-1)=B_n__2.value(11/2.0);
    A(n+1,n)=B_n__1.value(11/2.0);
    A(n+1,n+1)=B_n.value(11/2.0);

    b[0]=f_1;
    b[n+1]=f_n;

    Eigen::VectorXd a=Eigen::VectorXd::Zero(n+2);
    a=A.colPivHouseholderQr().solve(b);
    //std::cout<<"A："<<endl<<A<<endl;
    // 在求解完线性方程组后，将系数存入
    coefficients.resize(n+2);
    for (int i = 0; i <n+1; i++) {
        coefficients[i]=a[i];
       // std::cout<<"系数："<<coefficients[i]<<endl;
    }
}

double QuadraticBSplines::B_value(double x_value){
        double result=0.0;

        int n=knots.size();
        for(int i=0;i<n+2;i++){
            B_Basis_function B(2,i-1,knots);
          //  std::cout<<"B_value: "<<B.value(x_value)<<std::endl;
            result += coefficients[i]*B.value(x_value);
           // std::cout<<"coefficients: "<<coefficients[i]<<std::endl;
            
        }
       // std::cout<<"result: "<<result<<std::endl;
        return result;
}

#endif